8.4 s

Deriving a SIR epidemiological model

In this notebook we will go from a discrete epidemiological model to its continuous counterpart. We will go from a (discrete) system of difference equiations to a system of (continuous) differential equations. But first: SIR.

Susceptible, Infected, Recovered (SIR)

In this workhorse of epidemiology we model agents as belonging to a certain class, or compartment. Either they are

  • susceptible to get infected,

  • infected, or

  • removed/recovered

The models you may have heard about during the COVID19 pandemic are more sophisticated versions of what we will see here, but they share some important common features. At the end of this notebook, we will want to generate the following plot. It shows the number of agents belonging to each class at a certain point in time after the outbreak of the disease:

4.0 ms
37.0 ns

We will make this very simple. Instead of trying to model in great detail which type of agent (old, young, where, what time of the day, doing which activity etc) is likely to get infected, we will try to condense the overall probability of infection into a single number. Let's start to model the recovery from an infection IR:

we have N people infected at time 0. If each has probability p to recover each day, how many are still infected at day number t?

For any given person in I state, the probability to recover in a certain time step (i.e. each day) is a

3.2 μs
bernoulli (generic function with 1 method)
21.6 μs
N
200
27.0 ns
step! (generic function with 1 method)
28.1 μs
T
100
40.0 ns
83.3 μs
1
23.8 ms
5.8 s

ok great, now we know how many are I or R after t days with a certain probability p. Let's just condense this plot above and count how many infected there are each day!

8.3 μs
simulate_recovery (generic function with 1 method)
31.7 μs
652 ms
25.6 ms
515 ms
279 ms
372 ms

Deterministic dynamics for the mean: Intuitive derivation

20.5 μs

The mean seems to behave in a rather predictable way over time. Can we derive this?

Let It be the number of infectious people at time t. This decreases because some people recover. Since people recover with probability p, the number of people that recover at time t is, on average, pIt. [Note that one time unit corresponds to one sweep of the simulation.]

So

It+1=ItpIt

or

It+1=It(1p).

8.5 μs

At time t there are It infectious. How many decay? Each decays with probability p, so on average pIt recover, so are removed from the number of infectious, giving the change

ΔIt=It+1It=pIt

8.9 μs

We can rearrange and solve the recurrence relation:

It+1=(1p)It. so It+1=(1p)(1p)It1=(1p)2It1

and hence solve the recurrence relation:

It=(1p)tI0.

16.8 μs

Let's compare the exact and numerical results:

8.6 μs
exact
5.4 μs
65.1 ms

They agree well, as they should. The agreement is expected to be better (i.e. the fluctuations smaller) for a larger population.

8.7 μs

Continuous time

If we look at the graph of the mean as a function of time, it seems to follow a smooth curve. Indeed it makes sense to ask not only how many people have recovered each day, but to aim for finer granularity.

Suppose we instead increment time in steps of δt; the above analysis was for δt=1.

Then we will need to adjust the probability of recovery in each time step. It turns out that to make sense in the limit δt0, we need to choose the probability p(δt) to recover in time t to be proportional to δt:

p(δt)λδt,

where λ is the recovery rate. Note that a rate is a probability per unit time.

We get

17.2 μs

I(t+δt)I(t)λδtI(t)

2.5 μs

Dividing by δt gives

I(t+δt)I(t)δtλI(t)

We recognise the left-hand side as the definition of the derivative when δt0. Taking that limit finally gives

9.7 μs

dI(t)dt=λI(t)

That is, we obtain an ordinary differential equation that gives the solution implicitly. We know this has a general solution of the form

I(t)=Cexp(λt),

and together with an initial value of I(0)=I0 we get (just set t=0 and use the given value)

11.6 μs

I(t)=I0exp(λt).

Starting at time zero (t=0) and at a number of I0 infected, this is exponential decay at rate λ. Of course with a positive number (i.e. no minus sign, as here), we get exponential growth.

9.8 μs

SIR model

3.1 μs

Now let's extend the procedure to the full SIR model, SIR. Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with a certain probability.

3.1 μs

Let's denote by St and It be the number of susceptible and infectious people at time t, respectively, and by N the total number of people.

On average, in each sweep each infectious individual has the chance to interact with one other individual. That individual is chosen uniformly at random from the total population of size N. But a new infection occurs only if that chosen individual is susceptible, which happens with probability St/N. Then, upon meeting a susceptible person, the infection is transmitted only with a certain probability, b say.

Hence the change in the number of infectious people after that step is:

15.2 μs

ΔIt=It+1It=bIt(StN)

3.1 μs

It is useful to normalize by N, so we define

st:=StN;it:=ItN;rt:=RtN

4.3 μs

Including recovery with probability c we obtain the discrete-time SIR model:

5.1 μs

st+1=stbstitit+1=it+bstitcitrt+1=rt+cit

1.6 μs

And again we can allow the processes to occur in steps of length δt and take the limit δt0. With rates β and γ we obtain the standard (continuous-time) SIR model:

3.6 μs

ds(t)dt=βs(t)i(t)di(t)dt=+βs(t)i(t)γi(t)dr(t)dt=+γi(t)

1.7 μs

Note that no analytical solutions of these (simple) nonlinear ODEs are known as a function of time!

2.5 μs

Below is a simulation of the discrete-time model. Note that the simplest numerical method to solve (approximately) the system of ODEs, the Euler method, basically reduces to solving the discrete-time model! A whole suite of more advanced ODE solvers is provided in the Julia DiffEq ecosystem. We will in another notebook introduce that package, together with a continuous version of the SIR model.

7.8 μs
579 ms
0
38.0 ns
31.0 ns
34.0 ns
discrete_SIR (generic function with 2 methods)
46.3 μs
SIR
11.5 ms